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Abstract 

Fuel  cell  utilization  for  automobile  and  residential  applications  is  a  promising  option  in  order  to  help  reduce  environmental  concerns  such  as 
pollution.  However,  fuel  cell  development  requires  addressing  their  dynamic  behavior  to  improve  their  performances  and  their  life  cycle.  Since  the 
temperature  distribution  in  the  cell  is  known  to  be  an  important  factor  to  the  fuel  cell’s  efficiency,  a  cooling  device  is  often  added  to  homogenize 
the  temperature  in  the  cell  and  to  ensure  temperature  control. 

A  3D  dynamic  thermal  model  of  a  single  fuel  cell  is  presented  in  this  work  in  order  to  study  the  temperature  distribution  in  a  fuel  cell  cooled 
from  the  bottom  to  the  top  with  air.  The  model  is  governed  by  the  thermal  energy  balance,  taking  into  account  the  inlet  gas  humidity.  The  model  is 
developed  with  the  finite  difference  method  and  is  implemented  in  the  Matlab/Simulink  environment.  The  validation  is  based  on  the  performances 
of  the  “NEXA”  fuel  cell  produced  by  Ballard  Power  Systems. 

The  efficiency  analysis  of  that  air  cooling  device  reveals  that  the  cell  temperature  is  directly  linked  to  the  current  density  and  to  the  gas 
humidity — varying  from  30  °C  at  5 A  to  80  °C  at  35 A  at  low  humidity.  Moreover,  the  temperature  non-uniformity  in  the  stack  is  shown  to  be 
very  high.  As  a  result,  temperatures  are  higher  at  the  top  part  of  the  cell  than  at  the  bottom  part,  with  a  difference  of  up  to  a  5  °C.  Moreover  the 
non-uniformity  of  the  air  cooling  between  the  cells  of  the  stack  leads  to  large  temperature  variations,  up  to  8  °C,  from  one  cell  to  another.  These 
temperature  variations  result  in  large  voltage  disparities  between  the  cells,  which  reduce  the  total  electrical  power  of  the  entire  stack. 

©  2007  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

As  an  alternative  to  energy  resource  depletion  and  pollution 
problems,  fuel  cell  systems  have  received  much  attention  in 
research  and  development  during  the  last  few  years.  This  interest 
has  led  to  a  good  understanding  of  electrochemical  reactions  in 
fuel  cells.  Temperature,  which  is  an  important  factor  in  the  fuel 
cell’s  efficiency,  has  also  been  diversely  analyzed  and  modeled 
with  one  to  three-dimensional  thermal  models.  Since  tempera¬ 
ture  in  the  cell  influences  cell  humidity,  it  therefore  has  also  an 
indirect  influence  on  the  cell  output  power.  It  may  also  influence 
the  fuel  cell’s  durability  with  successive  drying  and  flooding.  A 
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detailed  understanding  of  the  stack  thermal  behavior  is  therefore 
necessary  for  the  choice  of  accurate  cooling  devices. 

As  noted  by  Berning  et  al.  [1],  the  temperature  distribution 
inside  the  fuel  cell  has  important  effects  on  nearly  all  trans¬ 
port  phenomena.  The  authors  emphasized  the  need  to  account 
for  thermal  gradients  and  multi-dimensional  transport  in  fuel 
cells  modeling.  The  effects  on  transport  have  also  been  analyzed 
by  Yan  et  al.  [2].  Moreover,  the  knowledge  of  the  temperature 
increase  due  to  irreversibilities  might  help  in  preventing  failures. 
As  reviewed  by  Biyikoglu  [3],  several  works  on  modeling  and 
analysis  of  fuel  cells  thermal  behavior  have  been  presented  in 
the  literature. 

Many  steady  state  models  have  been  published  [4-7].  Some 
[8,9]  assumed  a  uniform  temperature  throughout  the  cell  and 
others  [10-12]  present non-isothermal  models.  Meng  [10]  shows 
that  under  non-isothermal  two-phase  conditions,  the  heat  trans- 


K.P.  Adzakpa  et  al.  /  Journal  of  Power  Sources  179  (2008)  164-176 


165 


Nomenclature 

Am  membrane  effective  area  (m2) 

C  molar  concentration  (mol  m-3) 

Cp  specific  heat  capacity  (J  kg- 1  K- 1 ) 

£>A-b  diffusion  coefficient  of  the  gas  pair  A-B  (m2  s  - 1 ) 
Eq  thermodynamic  potential  at  standard  conditions 

(V) 

F  Faraday’s  constant  (C  mol- 1 ) 

h ,  /zfor,  hnSLt  convection  coefficients  (W  m-2  K-1) 

I  cell  load  current  (A) 

K,  Kx,  Ky ,  Kz  thermal  conductivity  coefficients 
(Wm-1K-1) 

n  number  of  transferred  electrons  in  the  reaction 

(two  for  H2  and  four  for  O2) 

N  molar  flow  (mol  m  2  s- 1 ) 
p  partial  pressures  (atm) 

Q  heat  source  (W) 

R  universal  gas  constant  (J  mol-1  K-1) 

RH  relative  humidity  (%) 

T  temperature  (K) 

V  voltage  or  overvoltage  (V) 

Greek  letters 
s  porosity 

k  water  content 

p  density  (kg  m- 3 ) 

o  conductivity  (Sm-1) 

Superscripts  and  subscripts 

act  activation 

ano  anodic 

cat  cathodic 

cel  relating  to  the  cell 

channel  relating  to  the  gas  channel 

CO  relating  to  CO 

GDL  relating  to  the  gas  diffusion  layer 

H2  relating  to  hydrogen 

in  inlet  quantity 

N2  relating  to  nitrogen 

O2  relating  to  oxygen 

ohm  ohmic  quantity 

out  outlet  quantity 

op  operational  quantity 

r  relating  to  a  reactant 

w  relating  to  water 


fer  process  significantly  increases  the  transient  response  time. 
In  his  model,  convection  effects  were  neglected.  Wu  et  al.  [12] 
presented  the  same  kind  of  results  which  were  highly  dependent 
on  the  membrane  thickness.  Nguyen  et  al.  [11]  also  addressed 
non-isothermal  behaviors  of  the  cell  from  the  point  of  view 
of  water  management  requirements  (humidification  or  water 
removal)  to  prevent  potential  membrane  dehydration  or  elec¬ 
trode  flooding.  They  showed  that  these  requirements  are  much 


more  conservative  than  those  predicted  assuming  isothermal 
conditions.  The  thermal  model  presented  by  Shan  and  Choe  [13] 
took  into  account  convection  as  well  as  conduction  phenom¬ 
ena.  The  authors  presented  a  one-dimensional  dynamic  thermal 
model  of  a  single  cell  based  on  a  mass  transport  description. 
They  analyzed  in  particular  the  effects  of  a  transient  tempera¬ 
ture  distribution  in  the  cell  on  its  electrical  performances.  Xue 
et  al.  [14]  also  presented  a  one-dimensional  dynamic  thermal 
model  in  which  the  convection  influences  are  the  most  impor¬ 
tant  factor.  In  their  model,  the  PEM  fuel  cell  is  subdivided  into 
three  control  volumes  representing  the  two  electrodes  and  the 
membrane.  Uniform  thermo-physical  properties  are  considered 
throughout  their  work. 

There  exist  several  CFD  based  fuel  cell  models  such  as 
[15, 16]  for  a  complete  analysis  of  thermal  transfer  in  a  single  cell 
of  PEMFC.  This  type  of  modeling  gives  a  detailed  microscopic 
description  of  the  whole  transport  mechanisms  in  the  cell.  The 
thermal  transient  response  is  analyzed  in  Ref.  [17].  Hwang  et 
al.  [18]  also  proposed  a  complete  PEMFC  thermal  model  based 
on  a  two-phase  temperatures  approach.  They  highlight  the  local 
thermal  non-equilibrium  between  the  solid  matrices  and  the  flu¬ 
ids.  Compared  to  all  these  models,  the  PEM  fuel  cell  model 
presented  in  this  paper  is  developed  for  the  purpose  of  control 
in  a  distributed  power  generation  system.  As  such,  the  model 
needs  to  be  fast  running;  it  ought  to  give  necessary  details  on  the 
thermal  behavior,  but  also  be  accurate  for  control. 

All  these  above  approaches  address  single  cell  modeling, 
although  cells  are  actually  assembled  in  stacks.  The  overall  heat 
balance  models  proposed  by  Koh  et  al.  [19]  and  Sohn  et  al.  [20] 
highlight  the  average  cell  temperature  variation  with  the  operat¬ 
ing  current  for  air  cooled  stacks.  They  showed  that  the  natural 
convection  is  not  sufficient  to  keep  the  cell  at  a  constant  temper¬ 
ature.  This  can  lead  to  large  cell  temperature  variations  resulting 
in  stack  voltage  degradation.  Additional  effects  of  local  temper¬ 
atures  on  the  voltage  and  power  losses  have  been  investigated 
by  Lee  et  al.  in  Ref.  [21]  with  a  dynamic  fuel  cell  stack  model. 
Other  dynamic  models  have  also  been  published.  For  instance, 
the  parametric  model  presented  by  Amphlett  et  al.  [22]  allows 
analyzing  the  dynamic  thermal  behavior  of  the  stack,  but  they 
assume  that  the  temperature  is  uniform.  Shan  and  Choe  [23] 
also  addressed  the  transient  thermal  distribution  in  the  stack. 
Their  models  characterize  the  dynamic  response  of  the  stack  in 
response  to  external  perturbations  (load  variations,  start-up  or 
shut  down,  etc.). 

The  main  goal  of  the  present  study  is  to  analyze  the  effect  of 
the  air  cooling  device  on  a  PEM  fuel  cell  stack  (NEXA  of  Bal¬ 
lard).  The  performances  of  this  stack  were  characterized  by  Zhu 
et  al.  [24]  and  del  Real  et  al.  [25]  by  experimental  measurements 
of  the  cell  voltage  as  a  function  of  gas  flow.  In  their  studies,  the 
effect  of  temperature  on  the  cell  voltage  is  emphasized,  although 
the  stack  temperature  is  assumed  uniform.  Our  experimental 
measurements  of  the  air  cooling  velocities  demonstrate  cooling 
non-uniformities  in  the  stack,  which  result  in  high  tempera¬ 
ture  disparities.  The  measured  air  velocities  help  to  define  the 
boundary  conditions  of  the  thermal  model.  A  three-dimensional 
dynamic  model  of  a  single  cell  is  developed  in  order  to  explain 
further  phenomena  such  as  the  cell  drying  or  flooding,  voltage 
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Fig.  1.  Description  of  the  modeled  PEM  fuel  cell  stack  and  its  ventilation. 


degradations,  etc.  The  developed  thermal  model  is  based  on  the 
total  thermal  power  produced  by  the  electrochemical  reactions. 
Whereas  the  thermal  model  is  three-dimensional,  the  electro¬ 
chemical  model  is  pseudo  two-dimensional;  it  helps  to  estimate 
heat  sources  in  the  cell.  The  heat  transfer  model  includes  the  con¬ 
duction  and  the  heat  generation  phenomena  inside  the  fuel  cell, 
and  the  convection  on  its  surface.  This  model  helps  to  understand 
how  the  natural  and  forced  convection  phenomena  influence  the 
temperature  field  in  the  cell.  The  thermal  results  are  presented 
and  analyzed  to  prove  the  influence  of  air  cooling  device  on  the 
fuel  cell  performances. 


2.  Model  description 

2.7.  Fuel  cell  description 


where  Gheat  (W),  Qah  (W),  Ge iec  (W)  and  Gsens  (W)  are 
respectively  the  thermal  power  to  be  removed  from  the  cell, 
the  theoretical  power  produced  by  electrochemical  reactions  of 
the  gases,  the  electrical  power  produced  by  the  fuel  cell,  and  the 
sensible  heat  of  the  gases.  The  last  term  2 sens  takes  into  account 
the  power  needed  to  heat  the  gases  to  the  cell  temperature. 

The  theoretical  power  supplied  to  the  cell  in  Eq.  (1)  by  the 
reacting  gases  is  given  by 


Qah  =  - 


A  77reac 
2  F 


(2) 


where  I  is  the  load  current  and  A77reac  is  the  total  enthalpy 
supplied  by  these  gases. 

The  electrical  power  produced  in  the  cell  is 


Gelec  —  kcell  ^ 


(3) 


The  aim  of  this  study  is  to  analyze  temperature  distribution 
in  typical  air  cooling  fuel  cells  and  its  impact  in  the  stack.  This 
analysis  is  based  on  the  “NEXA”  fuel  cell  produced  by  Ballard 
Power  Systems,  and  which  contains  47  vertical  cells  assembled 
in  series,  and  completely  separated  from  one  another  by  cooling 
channels  (Fig.  1).  The  stack  is  ventilated  from  the  bottom  to  the 
top.  Both  horizontal  surfaces  (tops  and  bottoms  of  the  cells)  are 
submitted  to  natural  convection  whereas  the  vertical  surfaces  are 
submitted  to  forced  convections. 

The  reference  axes  are  defined  as  follows  (Fig.  2): 

•  the  v-axis  is  along  the  cell’s  horizontal  width, 

•  the  y-axis  is  along  the  cell’s  vertical  depth  (cooling  channels 
length), 

•  and  the  z-axis  is  along  the  cell’s  thickness. 

2.2.  Global  heat  source  estimation 

In  a  running  fuel  cell,  a  large  part  of  the  energy  generated  is 
dissipated  as  heat  (approximately  50%  in  nominal  conditions). 
This  thermal  power  has  to  be  removed  from  the  cell  by  the 
cooling  device  to  avoid  important  overheating. 

The  thermal  power  dissipated  is  estimated  through  the  ther¬ 
modynamic  energy  balance  in  the  cell  [26] : 


Gheat  —  Qah  Gelec  H”  G sens 


(1) 


where  the  cell  voltage  Vceii  is  estimated  through  the  electro¬ 
chemical  model  presented  later  in  this  paper. 

The  sensible  heat  Gsens  represents  the  thermal  energy  vari¬ 
ation  of  the  gases  heated  in  the  cell.  It  is  estimated  from  the 
difference  between  the  thermal  energy  of  the  gases  entering  into 
the  fuel  cell  and  the  thermal  energy  of  the  unused  reaction  gases 
and  other  gases  leaving  the  fuel  cell  [22]. 

Gsens  =  Gs’ns  ~  Gsens  (4) 


Fig.  2.  Reference  system  for  the  single  cell  modeled:  (a)  A  cell  in  the  stack;  (b) 
A  single  ventilated  cell. 
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In  this  equation,  and  g°e ns  are  obtained  through  the  ther¬ 
modynamic  energy  balance  [27,28]  for  species  at  the  cell  inlets 
and  outlets,  respectively. 

el*  =  c,  (cPekn2  +  cPel:no  +  cPelnano) 

+  C  (cP  CS2  +  cP  Gn2  +  cp  2w,cat)  (5) 

/SOUt  T-'OUt  ( /-jOllt  I  Ollt  I  /-<  /-)OUt  \ 

£2sens  —  7ano  (CP^H2  cP^CO  +  cP^w,ano) 

+  P7  (cP  Q°£  +  cp  2n“‘  +  cP  Gwutcat)  (6) 

In  Eqs.  (5)  and  (6),  T™0,  7^t,  and  are  the  anode  and 
cathode  inlet  and  outlet  temperatures  (K).  <2J,nano,  0,  Q™ cat, 
2  s  cat  are  respectively  the  anode  and  cathode  inlet  and  outlet 
molar  fluxes  (mols-1)  of  species  s.  CVs  (Jmol-1  K_1)  is  the 
molar  specific  heat  capacity  of  species  s. 

2.3.  Dynamic  thermal  model 

Since  the  aim  of  the  dynamic  thermal  model  is  more  to  ana¬ 
lyze  the  influence  of  the  air  cooling  device  on  the  cell  than 
to  obtain  a  3D  thermal  description  of  the  cell,  some  simplify¬ 
ing  assumptions  have  to  be  made.  This  model  accounts  for  the 
dissipation  of  the  thermal  power  Qheat  from  the  cell  (see  Eq. 
(1)).  Within  the  cells,  only  conductive  heat  transfer  will  be  con¬ 
sidered,  since  many  authors  show  that  convective  heat  transfer 
inside  the  cell  is  negligible  [26] .  The  heat  removal  from  each  cell 
in  the  stack  is  ensured  by  forced  convection  on  the  lateral  sides 
(in  the  x-y  and  y-z  planes  of  Fig.  2)  and  by  natural  convection 
on  both  horizontal  surfaces  (top  and  bottom  in  the  x-z  plane). 
Because  of  the  relatively  low  temperatures  of  operation,  the  heat 
transfer  by  radiation  is  neglected. 

The  structure  of  the  stack  is  such  as  to  permit  the  assumption 
that  cell  temperatures  are  independent  from  one  cell  to  another. 
Indeed,  adjacent  cells  are  fully  separated  by  the  air  cooling  chan¬ 
nels.  In  this  case,  any  single  cell  of  the  stack  can  be  simulated 
with  its  corresponding  boundary  conditions  related  to  its  cooling 
channel.  Therefore,  the  whole  stack  is  described  as  the  assembly 
of  47  thermally  independent  cells. 

Conductive  heat  transfer  is  described  according  to  the  Fourier 
equation: 

v  •  (K  •  NT)  =  2 heat  (7) 

The  thermal  conductivity  coefficients  are  Kx,  Ky  and  Kz  in  the 
v,  y  and  z  directions,  respectively. 

Heat  transfer  between  the  cell  and  the  cooling  channels 
results  from  convective  exchanges.  The  boundary  condition  is 
expressed  as 

—K  •  VT  =  Ct"for  (8) 


The  thermal  energy  transmitted  by  forced  convection  to  the 
air  through  the  blades  is 

0foT  =  77/2 for (Tcell  -  Tair)  (9) 


where  27  is  the  blades’  efficiency,  hfOY  the  convection  coefficient 
and  Tair  is  the  air  temperature.  Tajr  is  estimated  from  experimen¬ 
tal  measurements  presented  in  the  following  section. 

The  convection  coefficient  is  estimated  with  the  following 
correlation: 


Nu 


Dhfor 

K 


(  Df1 3 

1.86  (  RePr—  J 


0.14 


(10) 


where  Nu,  Re  and  Pr  are  the  Nusselt  number,  the  Reynolds  num¬ 
ber  and  the  Prandtl  number,  respectively;  K and  D  are  the  thermal 
conductivity  coefficient  and  the  hydraulic  radius. 

With  air  (K^  =  0.03 Wm-1  K-1)  flowing  at  3 ms-1  in 
a  3.66  x  5.23  mm2  channel,  the  Nusselt  correlation  gives 
Nu  ex  5 . 1 ,  thus  /ifor  =  3 1 .3  W  m-2  K- 1 .  Thus,  the  forced  convec¬ 
tion  coefficient  is  set  to  hfOY  =  30  W  m-2  K-1  for  the  simulations. 

The  natural  convective  heat  flux  is  expressed  as  follows: 


0nat V  —  ^nat(7cell  T^mb)  (H) 

with  hnat  set  to  20%  of  hfor,  hnSLt  =  6  W  m-2  K_1 . 

These  values  are  in  good  agreement  with  the  mean  values  of 
the  natural  convection  coefficient  estimated  by  Koh  et  al.  [18] 
on  a  self-heated  PEMFC. 

In  this  model,  the  mathematical  cell  is  subdivided  into  finite 
volumes,  and  the  thermal  balance  in  each  control  volume  leads  to 
the  following  discreet  expression  for  an  internal  control  volume 
(i,j,k)  at  time  p 


rpp+ 1  _ 

iJ’k  ~  pCp 


+  ky 


+  k2 


Dx 2 

(rrP  _  OT'P  1  tP  1 

Vj-u  ZIij,k  +  1ij+ \,k) 

D p 

(rpP  OTP  _1_  tP  l 

LU  Jk-l  Z1iJ,k  +  Hj'k+l) 

D?  + 


(12) 


+  r 


i,j,k 


where  Dx,  D y,  D z  and  are  respectively  the  finite  element 
dimensions  along  the  x,  y  and  z  axes  and  in  time;  Cp  the  specific 
heat  capacity,  p  the  cell  density  and  q  is  the  thermal  power  pro¬ 
duction  rate.  Depending  on  the  position  of  the  control  volume 
(internal,  natural  convection  side  or  forced  convection  side),  27 
sub-cases  of  Eq.  (12)  are  defined. 

For  sake  of  thermal  stability  due  to  the  second  law  of  ther¬ 
modynamic,  Dt  is  bounded 


D  t  < 


_ pC  p _ 

2((^/Dx2)  +  ( KfDyb  +  (Kz/T>z2)  +  (hfm/Dx)  +  fe/D y)  +  (A^/Dz)) 


(13) 


where  p;ir  is  the  heat  flux  released  by  natural  or  forced  con-  where  hnat  and  hf0T  are  the  natural  and  forced  convection  coeffi- 
vection,  depending  on  the  considered  surface.  cients,  respectively. 
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2.4.  Cell  electrochemical  model 

As  presented  in  the  previous  section,  the  heat  released  in  a 
running  fuel  cell  is  directly  related  to  the  cell  voltage.  It  is  thus 
necessary  to  properly  compute  the  cell  voltage  to  estimate  the 
global  heat  source. 

Due  to  the  low  current  densities  in  the  modeled  stack,  it  is 
assumed  that  the  concentration  overvoltage  is  negligible,  and 
thus  the  cell  output  voltage  Vceii  is  given  by 


^cell  —  ^thermo  ^ohm  ^act  (14) 

The  thermodynamic  potential  Thermo  resulting  from  the  chem¬ 
ical  reaction  is  obtained  by  the  Nernst  equation  [26,29,30]: 

^thermo  =  Eq  (Top  —  A)) 


Fig.  3.  .  Description  of  the  modeled  PEM  fuel  cell. 


+  — f^ln  [(po2)l/2(Pti2)] 
nF 


(15)  2.5.  Mass  transport  model 


where  A S°  is  the  reaction  entropy  change  at  the  standard  tem¬ 
perature. 

The  ohmic  overvoltage  is 


l^ohm  —  ^ohm  1 


(16) 


where  the  total  resistance  in  the  cell  R0hm  depends  on  the  mem¬ 
brane  characteristics,  on  the  temperature  and  on  its  hydration 
levels.  It  is  given  by  the  following  integration  (17)  along  the 
membrane  thickness 

fLm  dx 

Rohm=  /  - - —  (17) 

Jo  AmcrTo  PW 

In  this  integration,  apop  (Sm-1)  is  the  membrane  conductiv¬ 
ity  at  the  operational  temperature  Top,  Am  (m2)  the  membrane 
effective  area  and  Lm  (m)  is  the  membrane  thickness.  The  mem¬ 
brane  conductivity  <jtov  depends  on  its  temperature  and  on  its 
water  content,  X,  expressed  as  the  molar  ratio  H2O/SO3-.  The 
membrane  conductivity  is  obtained  through  Eq.  (18)  [26,31]: 


7.6138  1.9796 

<%,  =  0.25  exp  (  4.1932  - 


-  1892  (  —  _ 

Top 


1 


353.15 


(18) 


The  water  content  distribution  along  the  membrane  thickness 
is  estimated  with  the  mass  transport  model  presented  below.  As 
a  result,  a  dryer  membrane  leads  to  higher  ohmic  losses. 

The  activation  overvoltage  is  defined  by  Eq.  (19)  with  a  semi- 
linear  parametric  transformation  of  the  Butler- Volmer  equation 
[32]: 

^act  =  §1  §2 Top  +  §3 Top  log  (cq2)  +  §4 Top  log(T)  (19) 

where  §1 ,  §2,  §3  and  §4  depend  on  the  cell  characteristics  [22,32]. 

The  electrochemical  model  is  strongly  dependent  on  the 
temperature  and  on  the  mass  transfer.  In  the  above  voltage  com¬ 
putation,  effective  reactants  partial  pressures  at  the  electrodes 
are  also  estimated  through  the  mass  transport  model. 


The  modeled  single  fuel  cell  is  fed  with  air  and  hydrogen 
and  is  schematically  represented  in  Fig.  3.  The  electrochemical 
model  is  a  pseudo  2D  model  in  which  the  anode  and  cathode 
channels  are  subdivided  into  control  volumes  where  the  reac¬ 
tant  concentrations  are  computed  in  order  to  obtain  their  mean 
concentration  at  the  catalyst  layers  (see  Fig.  3)  [26].  Since  the 
transient  responses  of  the  gases  are  clearly  lower  than  that  of 
liquid  water  in  the  membrane  (Dnq/DgSLS  a  10-4),  gas  transport 
is  assumed  to  be  in  steady  state.  Therefore,  the  continuity  equa¬ 
tion  in  steady  state  implies  that  all  the  mass  flows  are  uniform 
in  the  cell,  for  any  species  i 

VNi  =  0  (20) 


The  relationship  between  the  reactive  gas  flows  in  the  GDL  and 
the  load  current  are  given  by  Eq.  (21): 


(21) 


And,  because  the  membrane  is  impervious  to  gases,  the  inert  gas 
flow  (impurity  carbon  monoxide  at  the  anode  and  nitrogen  at  the 
cathode)  in  the  GDL  is  null.  The  vapor  water  flow  in  the  GDL 
is  estimated  with  the  mass  transport  model  in  the  membrane. 


2.5.1.  Gas  transport  in  the  gas  diffusion  layer  and  the 
bipolar  plates 

Along  their  progression  in  the  feeding  channels  of  the  bipo¬ 
lar  plates,  reactive  gases  are  consumed  at  the  electrodes  and  the 
resulting  water  must  be  evacuated.  To  model  this  transverse  gas 
transport,  two  flow  regimes  are  considered:  the  convective  flow 
of  the  gases  from  the  center  of  the  channel  to  the  porous  gas  dif¬ 
fusion  layer  (GDL),  and  the  diffusion  flow  of  the  gases  through 
the  GDL  to  reach  the  catalyst  layer. 

The  masses  exchanged  between  the  GDLs  and  the  channels 
are  modeled  with  a  convective  mass  exchange  coefficient 

Ns  =  k(CfDL  -  cf annel)  (22) 

where  k  (ms-1)  is  the  mass  transport  coefficient,  and  CpDL 
and  C£hannel  (molm-3)  are  respectively  the  concentration  of 
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species  s  at  the  GDL/channel  interface  and  in  the  feeding 
channel. 

In  the  GDL,  only  diffusive  transport  is  considered  because  of 
the  low  gases  velocities.  The  diffusive  mass  flow  is  expressed  as 

dcs 

Ns  =  -Dw-s-^  (23) 

dz 

where  Dw_s  is  the  water  vapor  diffusion  coefficient  in  species 
(gas)  s. 

The  vapor  water  concentration  at  the  GDL/membrane  inter¬ 
face  depends  on  the  membrane  water  content  via  the  sorption 
isotherm.  The  diffusion  coefficients  are  computed  using  the  Slat¬ 
tery  and  Bird’s  approach  [33]  and  Bruggemann’s  correction  [34] 
in  order  to  take  into  account  the  GDL  porosity. 


where  Da~b  (m2s-1)  is  the  diffusion  coefficient  for  the  AB 
binary  gas  mixture,  p  (atm)  the  mixture  pressure,  s  the  porosity 
of  the  electrode  and  b ,  c  and  d  are  gas  pair-related  constants. 

The  average  concentration  is  computed  for  each  species  in 
the  channels  in  order  to  estimate  the  mean  reactants  partial 
pressures/concentrations  at  the  reaction  site.  Consequently,  the 
effective  reactants  partial  pressures  are  assumed  uniform  on  the 
membrane’s  surface. 

Finally,  the  partial  pressures  pr  of  reactants  r  are  computed 
from  their  concentration  in  Eq.  (25): 


Pr  —  CrRTQ 


op 


(25) 


2.5.2.  Water  transport  in  the  membrane 

The  cathodic  active  layer  is  assumed  to  be  part  of  the  mem¬ 
brane  and  have  a  measurable  thickness.  The  water  mass  balance 
gives 


acw  dNw 

dt  +  dz 


Rh2o  in  the  cathode 
0  elsewhere  in  the  membrane 


(26) 


where  Nw  is  the  water  molar  flow  (mol m  2  s  ]). 

/?h20  is  the  water  generation  rate  (molm-3  s-1)  given  by 


Kh2o  = 


2AmLcatF 


(27) 


where  Lcat  is  the  cathode  thickness  (m),  Am  the  cell  section  (m2) 
and  i  is  the  current  density  (Am-2). 

The  water  transport  in  the  membrane  is  caused  by  three  phe¬ 
nomena:  diffusion,  electro-osmotic  drag  and  convection.  Water 
molar  flow  is  thus  expressed  as  [1 1] 


Nw  —  Dw 


9Cw 

dz 


+  Od  T-  CWVy 
F 


(28) 


where  Dw  is  the  water  diffusion  coefficient  (m2s-1),  r]d  the 
electro-osmotic  drag  coefficient  and  vw  is  the  water  velocity 
(ms-1).  The  water  diffusion  coefficient  depends  on  the  operat¬ 
ing  temperature  and  on  the  membrane  hydration  level  [3,1 1]. 

The  water  drag  coefficient  expresses  the  number  of  water 
molecules  pulled  by  each  proton  from  the  anode  to  the  cathode. 


For  Nation®  membranes,  the  following  expressions  are  often 
used  [9] 


2.5 

rid  =  — X 
1  22 


EW 

X  —  - 

Pm 


(29) 

(30) 


where  X  is  the  membrane  water  content  (mol  H2 O/mol  SO3-), 
EW  is  the  membrane  equivalent  weight  (membrane  mass 
(kg)/mol  SO3-),  and  pm  is  the  density  of  the  dry  membrane 
(kgm-3).  The  membrane  swelling  is  neglected. 

By  neglecting  gravity  and  assuming  a  linear  pressure  drop 
across  the  membrane,  the  water  velocity  is  given  by  [11] 


KdP  =  K  /p&-pc\ 
M  dz  fi  \  Lm  ) 


(31) 


where  K  is  the  membrane  permeability  (m2),  pi  the  water  vis¬ 
cosity  (Pa  s),  pa  and  pc  are  respectively  the  absolute  pressure  on 
the  anode  and  cathode  side  of  the  membrane  (Pa)  and  Lm  is  the 
membrane  thickness  (m). 

The  following  second  order  partial  derivative  equation  is 
obtained  by  replacing  (28)— (31)  in  Eq.  (26): 


9Cw 

dt 


d 

dz 


dCw 
'  dz 


ndC% 

P^L  + 

dz 


0 

Ph2o 


with 


2.5  i  EW 
22  F  pm 


and 


(  Pa  ~Pc  \ 

V  Lm  ) 


(32) 


Eq.  (32)  is  then  solved  using  an  implicit  upwind  finite  difference 
scheme  [1 1]  to  obtain  the  water  concentration  Cw  and  the  water 
flows  Nw  in  the  whole  membrane  (via  Eq.  (28)).  The  bound¬ 
ary  conditions  are  the  water  concentration  C^no  and  C^at  at  the 
membrane  interfaces.  The  water  concentrations  are  converted 
into  water  content  using  Eq.  (30). 


2.5.3.  Mass  balance  along  the  channel 

For  numerical  simulation  of  the  reactant  gas  flow  (in  the  y- 
direction — see  Fig.  3),  the  feeding  channels  are  subdivided  into 
m  control  volumes  (40  for  the  oxygen  channel  at  the  cathode 
and  10  for  the  hydrogen  channel  at  the  anode)  in  which  the  mass 
balances  are  computed.  Reactants  inlet  conditions  are  applied  to 
the  first  control  volume  whereas  the  conditions  in  the  last  control 
volume  correspond  to  those  of  the  cell  outlet. 

The  concentration  Clr+l  of  reactant  r  in  the  (i  +  l)th  control 
volume  is 

c;+l  =  cl-—  -5L_  with =  (33) 

n  F  ^channel  m  Vg 

where  A  ns  the  time  required  for  the  gas  flow  to  get  across  the 
control  volume,  Vchannel  (m3)  the  total  volume  of  the  channel, 
/-channel  (m)  the  length  of  the  channel  and  vg  (ms-1)  is  the  gas 
velocity  along  the  channel. 
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Similarly,  the  vapor  water  concentration  is  computed  through 
the  following  equation: 


W+l  _  r'i 

— 


<+1Af 

^channel/  m 


(34) 


Reactants  mean  concentrations  are  then  computed  in  the  chan¬ 
nel.  The  current  density  is  assumed  to  be  uniform  throughout 
the  cell.  So,  given  the  mean  concentration  C^ean  (molm-3)  of 
reactant  r  in  the  channel,  its  effective  concentration  Cr  at  the 
reaction  site  can  be  computed  by  the  following  equation: 


_  £*mean  _  ^  2-Pw— x  ~l~ 

;  1  nF  Dw-rAmk 


(35) 


where  n  is  equal  to  two  and  four  for  hydrogen  and  oxygen 
respectively. 

The  flow  of  vapor  water  transferred  from  the  electrode  to  the 
gas  channel  in  the  (i  +  l)th  control  volume  is  computed  according 
to  the  following  equation: 


Ni+1  = 

1  Y  \\7 


Dw—xAftk 

2£>w-x  +  Tgdl^ 


(CJ 


i+t 


Ce ) 


(36) 


where  Z)w_ g  is  the  water  vapor  diffusion  coefficient  in  gas  g , 
Ad  the  membrane  area  encompassed  by  the  control  volume,  k 
(ms-1)  the  mass  transport  coefficient,  Lgdl  (m)  the  GDL  thick¬ 
ness,  and  (molm-3)  is  the  water  vapor  concentration  at  the 
electrode. 


3.  Results  and  discussion 

3.1.  Experimental  set-up 

In  the  PEM  fuel  cell  thermal  model,  each  cell  has  the  geom¬ 
etry  shown  in  Fig.  4. 

The  grey  part  is  composed  of  a  0.2  mm  thick  membrane,  two 
0.05  mm  thick  catalyst  layers  and  two  0.5  mm  thick  diffusion 
layers.  The  equivalent  thermal  properties  of  the  whole  cell  are 
estimated  with  the  thermal  properties  of  each  material,  according 
to  the  literature  review  presented  in  [35] .  Thus,  the  through-plane 
thermal  conductivity  of  the  cell  is  set  to  Kz  =  0.17  W  m-1  K-1. 
The  in-plane  thermal  conductivity  ( Kx  =  Ky  =  0.85  W  m-1  K-1) 


Catalyst  layers+ 
Diffiisinn  laver 

Fig.  4.  Single  cell  geometry. 


is  chosen  five  times  greater  than  the  through-plane  one’s  to 
represent  the  thermal  anisotropy  of  these  materials  [35,36]. 
Moreover,  the  cell  has  an  average  density  of  2200  kg  m- 3  and 
a  710  J  kg-1  K-1  specific  heat  capacity.  The  two  other  vertical 
parts  represent  the  1.6  mm  thick  blades  and  are  assumed  to  be 
in  aluminium  [7].  Their  density  and  specific  heat  capacity  are 
2707  kg  m-3  and  896  J  kg-1  K-1,  respectively. 

According  to  the  heat  source  estimation  presented  previously, 
the  thermal  power  to  remove  from  the  cell  versus  the  load  cur¬ 
rent  is  plotted  in  Fig.  5  with  10%  and  90%  relative  humidity  of 
the  inlet  gases.  In  Fig.  5,  RHano  and  RHcath  are  respectively  the 
anode  and  cathode  inlet  gas  relative  humidity.  For  any  current, 
the  sensible  heat  ( <2 sens)  is  less  than  1  W  in  both  cases  and  has 
a  reduced  influence  on  the  thermal  power  to  remove  from  the 
cell.  On  the  other  hand,  the  gas  humidification  has  a  significant 
impact  on  fuel  cell  performances  due  to  its  large  effect  on  the 
membrane  resistance.  The  thermal  power  versus  the  load  cur¬ 
rent  is  summarized  in  Fig.  5  for  low  and  high  relative  humidity 


ForRH  =10%  and  RH  =10%  ForRH  =90%  and  RH  =90% 

ano  cath  ano  cath 


Fig.  5.  Cell  thermal  power  vs.  load  current  rceii  =  70  °C. 
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(a)  Air  speed  with  high  ventilation 


Air  speed  (m/s) 


A  D 

(b)  Air  speed  with  medium  ventilation 


Fig.  6.  Experimental  air  flow  on  top  of  the  cell. 


At  high  relative  humidity  (RHano  =  RHcath  =  90%),  the  theoret¬ 
ical  power  splits  into  two  comparable  parts:  the  part  converted 
into  electricity  and  the  part  released  as  heat  (see  Eq.  (1)).  More¬ 
over,  at  high  currents,  the  energy  released  as  heat  (more  than 
27  W  at  35  A)  becomes  more  important  than  the  energy  con¬ 
verted  into  electricity  (around  23  W  at  35  A).  On  the  other  hand, 
at  low  gas  humidification  (RHano  =  RHcath  =  10%),  the  electrical 
power  is  restricted  to  14  W  at  35  A  because  of  the  high  mem¬ 
brane  resistance  (up  to  0.78  Q  cm2  compared  to  0.06  Q  cm2  with 
90%  relative  humidity  of  the  inlet  gases).  As  a  result,  the  power 
released  as  heat  in  the  cell  increases  up  to  38  W  at  35  A.  These 
figures  show  clearly  that,  with  the  exception  of  £9  sens,  the  other 
three  powers  not  only  increase  with  the  current,  but  that  Qheat 
and  2  elec  depend  noticeably  on  the  inlet  gas  humidity.  The  sim¬ 
ulation  of  the  temperature  profile  in  the  cell  is  based  on  this 
model  of  the  computed  thermal  power. 

3.2.  Boundary  conditions  of  the  cell  thermal  model 

In  order  to  analyze  the  efficiency  of  the  cooling  device,  the 
air  velocity  was  measured  at  the  top  of  the  cell  stack,  in  the 
plane  ABCD  in  Fig.  1  (where  the  stack  length  (AD,  BC)  is 
56  cm  and  the  stack  width  (AB,  DC)  is  25  cm)  in  maximum  and 
medium  ventilation  conditions  (Fig.  6 — the  numbers  along  the 


color  scale  on  the  right  represent  air  velocity).  The  air  velocity 
was  measured  with  a  mobile  anemometer  in  42  zones  uniformly 
distributed  over  the  top  of  the  stack.  The  top  surface  is  subdivided 
into  seven  regions  along  the  stack  length  and  six  along  the  stack 
width. 

Fig.  6  shows  that  the  air  flow  is  not  uniform  throughout  the 
stack;  some  channels  have  a  faster  air  flow  than  others.  A  com¬ 
parison  of  Fig.  6a  and  b  shows  that  the  zones  of  highest  and 
the  lowest  air  flow  remain  in  the  same  place  even  if  the  ventila¬ 
tion  changes.  The  boundary  conditions  are  therefore  chosen  to 
take  this  cooling  non-uniformity  into  account.  Different  cooling 
channels  are  analyzed  in  order  to  represent  the  thermal  behavior 
of  the  whole  stack. 

The  boundary  conditions  in  the  thermal  model  are  given  by 
the  temperatures  measured  in  the  cooling  channel  along  the  y- 
axis  (vertical  channel  length  axis).  In  order  to  analyze  the  thermal 
behavior  of  each  cell  in  the  stack,  five  sample  channels  have  been 
chosen  to  analyze  the  temperature  profile  of  the  whole  stack.  The 
positions  of  the  sample  channels  are  shown  in  Fig.  7a.  Channel 
1  represents  the  lowest-cooled  (due  to  lowest  air  flow)  channel, 
whereas  channel  5  represents  the  highest-cooled  channel.  For 
each  of  these  channels,  the  temperature  profiles  are  measured 
along  the  channel  length.  The  air  temperature  measurements  are 
used  as  boundary  conditions  in  the  thermal  model.  The  tem- 


Temperature  as  a  function  of  channel  depth  with  a  5A  load  current  Channel  1 


(a) 


(b) 


(c) 


Fig.  7.  Boundary  conditions  in  sample  channels:  (a)  Sample  channel  positions,  (b)  Channel  temperature  profile  (with  a  35  A  current)  and  (c)  Channel  temperature 
profile  (with  a  5  A  current). 
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perature  profiles  in  these  channels  for  35  A  and  5  A  currents,  as 
shown  respectively  in  Fig.  7b  and  c,  give  the  boundary  conditions 
for  the  numerical  analysis. 

The  measured  temperatures  in  the  channels  allow  a  polyno¬ 
mial  analytical  estimation  of  the  temperature  profiles  along  the 
channels  length  (Fig.  7b  and  c). 

As  expected,  the  temperature  increases  along  the  cooling 
channels  from  bottom  (right)  to  top  (left)  as  air  is  heated  inside 
the  cell.  The  difference  between  outlet  and  inlet  temperatures 
is  typically  about  3  °C  at  low  current  (Fig.  7b:  channel  1  tem¬ 
perature  increases  from  28.3  °C  to  31.3  °C),  while  it  is  about 
23  °C  at  high  current  (Fig.  7c:  channel  1  temperature  increases 
from  33  °C  to  55  °C).  At  high  currents,  the  air  flow  is  not  suffi¬ 
cient  for  a  convenient  cooling  of  the  cells,  even  though  the  fan 
is  internally  controlled  by  the  system  (higher  air  flow  at  higher 
load  current).  Such  a  large  temperature  gradient  will  result  in  a 
non-uniform  temperature  distribution  throughout  the  cell.  More¬ 
over,  at  this  high  current,  the  non-uniform  air  ventilation  leads  to 
a  5  °C  difference  in  air  temperature  between  the  channels,  pro¬ 
ducing  relatively  great  temperature  non-uniformities  throughout 
the  stack  (the  corresponding  figure  at  low  current  is  about  1  °C). 
The  use  of  water  instead  of  air  as  a  cooling  fluid  should  reduce 
the  temperature  difference  between  inlet  and  outlet,  since  the 
specific  heat  of  water  is  some  4000  times  greater  than  that  of 


air.  This  would  produce  smaller  temperature  gradients  along  the 
cells,  resulting  in  smaller  temperature  variations  between  cells. 

3.3.  Steady -state  thermal  fields 

Fig.  8  shows  the  computed  simulation  of  the  steady- state  tem¬ 
perature  field  of  a  cell  operated  with  a  35  A  load  current  and  at 
90%  relative  humidity,  and  submitted  to  the  lowest  cooling  inten¬ 
sity  (channel  1).  Fig.  8a-d  show  respectively  the  cell  external 
temperature,  the  temperature  field  in  the  middle  of  the  cell  width 
(in  the  y-z  plane),  the  temperature  field  in  the  middle  of  the  cell 
thickness  (x-y  plane)  and  the  temperature  field  in  the  middle  of 
the  cell  depth  (x-z  plane).  The  lowest  and  highest  temperatures 
in  the  cell  are  67  °C  and  70  °C  respectively.  The  results  show 
that  along  the  y-axis,  the  temperature  decreases  from  the  top  to 
the  bottom  because  the  cells  are  cooled  from  the  bottom  to  the 
top.  As  expected,  results  also  show  that  the  center  of  the  cell, 
including  the  membrane  and  the  catalyst  and  diffusion  layers,  is 
hotter  due  to  the  thermal  power  generation  inside  the  cell.  The 
temperature  field  is  symmetrical  along  the  cell  thickness  due  to 
the  symmetrical  cooling  conditions  on  both  sides  of  the  cell. 

It  should  be  emphasized  that  the  main  purpose  of  this  paper 
is  more  to  analyze  the  influence  of  the  cooling  device  on  the 
thermal  distribution  (from  bottom  to  top)  along  the  cell,  than 


Fig.  8.  Steady-state  cell  temperature  field  -  35  A  load;  current  -  channel  1 :  (a)  Cell  external  temperature,  (b)  A  cut  at  x  =  58  mm,  (c)  A  cut  at  z  =  2.2  mm  and  (d)  A 
cut  at  y  =  60  mm. 
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to  analyze  internal  3D  temperature  profiles  within  the  cell.  This 
will  help  to  anticipate  its  humidification  (flooding  and  drying).  In 
order  to  study  the  effect  of  air  cooling  magnitude,  Fig.  9  repeats 
the  simulation  of  Fig.  8  (cell  operating  at  35  A  and  at  90%  rel¬ 
ative  humidity),  but  this  time  submitted  to  the  highest  cooling 
intensity,  as  seen  in  channel  5.  In  this  case,  the  extreme  temper¬ 
atures  are  59  °C  and  62  °C  (compared  to  67.5  °C  and  70  °C  in 
the  low-magnitude  cooling  of  Fig.  8c).  It  thus  appears  that  an 
increase  of  the  cooling  intensity  does  not  have  a  strong  effect 
on  the  internal  cell  temperature  gradient.  However,  due  to  the 
better  cooling,  the  fuel  cell  temperature  decreases  by  about  8 
degrees.  Furthermore,  the  same  observations  about  temperature 
fields  and  symmetries  as  in  the  previous  case  can  be  done. 

By  comparing  the  cell  temperatures  presented  in  Figs.  8  and  9, 
it  appears  that  for  a  35  A  operating  current  and  a  90%  relative 
humidity,  a  temperature  variation  of  1 1  °C  throughout  the  stack 
can  be  expected:  3  °C  arising  from  the  internal  cell  gradient 
and  8  °C  from  the  temperature  difference  between  the  coolest 
and  the  warmest  cells.  This  means  that  the  different  cells  of  the 
stack  may  have  different  electrochemical  behavior  due  to  their 
different  temperatures.  This  remark  is  confirmed  by  [24],  who 
reveals  cell  voltage  variation  in  the  stack,  especially  for  the  two 
cells  nearest  to  the  air  compressor. 


In  the  same  electrical  operating  conditions  as  above  but 
with  a  10%  relative  humidity  only,  the  shape  of  the  tempera¬ 
ture  profile  remains  the  same.  However,  the  cell  temperature 
ranges  between  87  °C  and  92  °C  in  the  lowest  cooling  con¬ 
ditions  and  between  80  °C  and  84  °C  in  the  highest  cooling 
conditions.  As  expected,  the  cell  temperature  increases  due 
to  the  higher  membrane  resistance  resulting  in  higher  ohmic 
losses.  The  temperature  gradient  inside  the  cell  therefore  also 
increases  up  to  5  °C  at  10%  relative  humidity  instead  of  3  °C 
with  90%  relative  humidity.  Throughout  the  stack,  the  tem¬ 
perature  gradient  due  to  the  non-uniform  air  cooling  changes 
from  11°C  at  90%  relative  humidity  to  12  °C  at  10%  rela¬ 
tive  humidity.  Therefore,  the  influence  of  the  relative  humidity 
is  significant  in  the  cell,  but  is  negligible  from  one  cell  to 
another. 

The  steady- state  temperature  field  of  a  cell  operated  at  5  A 
was  also  studied.  These  results  show  that  the  predicted  temper¬ 
atures  of  the  cell  are  fairly  uniform  around  32  °C,  as  measured 
experimentally  by  [24].  The  temperature  gradient  does  not 
exceed  0.5  °C  throughout  the  cell.  It  can  therefore  be  concluded 
that  the  higher  the  operating  current,  the  higher  the  cell  temper¬ 
ature  and  the  more  important  the  temperature  gradients  within 
the  stack. 


Fig.  9.  Steady  state  cell  temperature  field  -  35  A  load  current  -  channel  5:  (a)  Cell  external  temperature,  (b)  A  cut  at  x  =  58  mm,  (c)  A  cut  at  z  =  2.2  mm  and  (d)  A  cut 
at  y  =  60  mm. 
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3.4.  Transient  state  thermal  profile 

The  transient  thermal  response  of  the  cell  submitted  to  a  load 
current  demand  variation  from  5  A  to  35  A  with  a  90%  rela¬ 
tive  humidity  is  plotted  in  Fig.  10,  according  to  the  boundary 
conditions  given  by  channel  1  (the  highest  temperature  zone). 

In  this  figure,  the  transient  temperature  in  the  centre  of  the  cell 
is  shown  in  Fig.  10a.  Transient  temperature  profiles  are  plotted 
along  three  significant  axes:  along  the  cell  thickness  z,  along  the 
cell  depth  y,  and  along  the  cell  width  v. 

The  same  conclusions  concerning  temperature  distributions 
can  be  drawn  as  for  the  steady  state  case.  Moreover,  along  the  cell 
thickness  (Fig.  10b),  the  temperature  difference  increases  during 
the  transient  phase  -  because  of  the  thermal  power  generation 
in  the  heart  of  the  cell  -  and  reaches  its  maximum  of  2  °C  at 
steady  state  (temperature  at  1500  s  in  Fig.  10b).  Along  the  cell 
depth  (Fig.  10c),  the  temperature  difference  also  increases  in 
the  transient  phase  and  reaches  its  maximum  of  3  °C  at  steady 
state  (temperature  at  1500  s  in  Fig.  lOd).  This  is  due  to  the 
ventilation  going  from  the  bottom  to  the  top  of  the  cell,  thus 
implying  a  better  cooling  of  the  bottom  than  the  top  of  the 
cell. 

Even  though  the  transient  phase  duration  is  up  to  1500  s,  it 
can  be  seen  from  Fig.  10a  that  the  temperature  reaches  99%  of 
its  steady  state  value  after  900  s.  These  transient  responses  are 
pretty  slow,  so  that  thermal  steady  state  should  not  be  reached 
easily  in  transient  operation  of  the  cell. 


Transient  temperature  (°C)  at  A  cut  atX=58mm.  Y=60mm  and  Z=2.2mm 


3.5.  Effect  of  the  temperature  on  voltage 

This  study  looks  at  the  temperature  variations  in  air  cooled 
stacks,  depending  on  the  operating  conditions.  At  low  currents, 
the  cell  temperature  is  low  (around  30  °C)  and  pretty  uniform 
on  the  cell  surface.  However,  an  increase  of  the  current  leads 
to  the  increase  of  the  mean  cell  temperature  (up  to  75  °C  with 
a  90%  relative  humidity  and  up  to  92  °C  with  a  10%  relative 
humidity,  at  35  A);  it  also  leads  to  the  increase  of  the  local  non¬ 
uniformity,  as  observed  experimentally  by  [37].  As  a  result  of 
this  large  temperature  range,  the  electrical  behaviors  of  the  cell 
change.  The  influence  of  the  local  temperature  variation  on  fuel 
cell  voltage  is  emphasized  in  several  publications  [12,38].  Actu¬ 
ally,  mass  and  charge  transfers  strongly  depend  on  temperature. 
The  results  also  show  that  the  gas  humidities  have  a  great  impact 
on  heat  sources,  and  on  temperature  disparities  in  the  cell.  Low 
gas  humidity  results  in  high  cell  temperatures.  And  as  the  vapor 
saturation  pressure  depends  on  temperature,  an  increase  in  tem¬ 
perature  reduces  the  gas  relative  humidity.  This  involves  dryness 
of  the  membrane  and  consequently  degradation  in  the  cell  volt¬ 
age.  All  together  this  means  that  running  a  fuel  cell  with  low 
inlet  gas  humidity  can  result  in  membrane  failure  if  the  cell  is 
not  properly  cooled,  especially  at  high  current  densities. 

On  the  other  hand,  high  gas  humidities  can  result  in  flooding 
of  the  active  layers,  GDL  and/or  feeding  channels.  In  the  above 
temperature  ranges,  the  saturation  pressure  is  pretty  low  (only 
0.043  bar  at  30  °C);  water  condensation  may  therefore  occur 


Fig.  10.  Cell  transient  temperature  profile  -  35  A  current  -  with  a  5  A  current  as  initial  condition:  (a)  Temperature  in  the  centre  of  the  cell,  (b)  A  cut  at  X=  58  mm 
and  Y=  60  mm,  (c)  A  cut  at  X = 58  mm  Z = 2.2  mm  and  (d)  A  cut  at  Y=  60  mm  and  Z  =  2.2  mm. 
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in  the  cell.  A  discussion  about  the  preferential  sites  of  water 
condensation  is  presented  in  [34]  and  descriptions  of  two-phase 
flows  in  a  non-isothermal  cell  can  be  found  in  [39].  The  authors 
insist  on  liquid  water  accumulation  in  the  cell  which  drastically 
reduces  the  cell  voltage  at  high  currents  and  high  gas  humidity. 
Thus,  the  simulated  cell  potentials  are  intentionally  not  plotted 
here  because  the  model  does  not  yet  take  into  account  the  water 
condensation  phenomenon.  Diphasic  transport  of  water  in  the 
cell  is  currently  under  study.  These  remarks  lead  to  the  following 
statement  for  the  air  cooling  device  studied  here:  at  high  currents, 
the  temperature  difference  between  the  top  and  the  bottom  of 
the  cell  could  lead  to  membrane  drying  at  the  top  because  of  the 
higher  temperatures,  and  to  water  condensation  at  the  bottom 
due  to  the  lower  temperatures.  Moreover,  water  condensation 
at  the  bottom  would  be  increased  due  to  the  gravity.  So,  large 
voltage  losses  can  be  encountered  in  the  whole  cell  because  of  the 
temperature  non-uniformity.  Moreover,  the  disparities  in  cooling 
air  velocities  from  one  cell  to  another  are  also  pointed  out.  This 
implies  that  cells  in  the  stack  have  different  electrochemical 
behaviors  depending  on  their  temperature.  Thus,  some  cells  in 
the  stack  are  more  subjected  to  drying  while,  at  the  same  time, 
some  others  are  subjected  to  flooding.  It  follows  that  the  cooling 
device  for  the  NEXA  stack  has  to  be  improved  to  keep  the  cell 
at  its  optimal  temperature. 

4.  Conclusions 

The  aim  of  this  study  is  to  characterize  the  influence  of  the  air 
cooling  device  on  fuel  cell  performances.  To  reach  this  goal,  a 
dynamic  3D  thermal  modeling  of  a  proton  exchange  membrane 
fuel  cell  is  presented.  The  model  includes  the  natural  convec¬ 
tion  due  to  ambient  air  cooling,  the  forced  convection  due  to  the 
stack  ventilation,  and  the  influence  of  the  thermal  power  gener¬ 
ation  inside  each  cell.  The  gas  humidity  shows  a  great  impact 
on  membrane  resistance,  resulting  in  higher  heat  generation  in 
the  cell.  Different  air  flow  zones  are  studied  to  understand  the 
thermal  behavior  of  the  cells.  This  model  brings  out  the  temper¬ 
ature  non-uniformity  throughout  the  stack.  The  “NEXA”  fuel 
cell  stack  of  Ballard  Power  Systems  was  used  to  validate  the 
approach. 

In  the  cell  depth  direction,  the  temperature  difference 
between  the  top  and  the  bottom  of  the  cell  reaches  its  maximum 
of  5  °C  at  steady  state  at  35  A.  This  difference  is  linked  to  the  air 
ventilation  from  the  bottom  to  the  top  of  each  cell.  Moreover,  the 
influence  of  the  load  current  on  the  mean  cell  temperature  is  also 
emphasized.  For  instance,  the  temperature  of  the  cell,  fed  with 
high-humidified  gases,  is  around  30  °C  for  a  5  A  load  current  and 
around  70°Cfora35A  load  current.  With  low  gas  humidity, 
the  total  heat  production  rises,  resulting  in  higher  temperatures 
in  the  cell  (up  to  92  °C  at  35  A).  The  results  also  demonstrate 
temperature  non-uniformity  in  the  stack,  which  increases  with 
the  load  current:  3-5  °C  difference  in  the  cell  and  up  to  7-8  °C 
difference  between  the  cells,  due  to  the  cooling  non-uniformity 
at  35  A.  Finally,  the  influence  of  such  temperature  differences  on 
electrical  performances  is  analyzed.  The  temperature  distribu¬ 
tion  is  responsible  for  electrical  disparities  in  the  cell  depth  and 
from  one  cell  to  another.  Moreover,  the  temperature  field  tends  to 


involve  water  condensation  at  the  bottom  and  membrane  dehy¬ 
dration  at  the  top.  This  results  in  cell  voltage  disparities,  which 
reduces  the  global  electrical  power  produced  by  the  stack.  This 
characterization  will  be  useful  for  an  efficient  control  of  the  fuel 
cell  stack. 
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